#include "aimage.h"
#include "aimage_dcraw.h"
#include "astar_file.h"
#include "aimage_ppm.h"

/*
   dcraw.c -- Dave Coffin's raw photo decoder
   Copyright 1997-2010 by Dave Coffin, dcoffin a cybercom o net

   This is a command-line ANSI C program to convert raw photos from
   any digital camera on any computer running any operating system.

   No license is required to download and use dcraw.c.  However,
   to lawfully redistribute dcraw, you must either (a) offer, at
   no extra charge, full source code* for all executable files
   containing RESTRICTED functions, (b) distribute this code under
   the GPL Version 2 or later, (c) remove all RESTRICTED functions,
   re-implement them, or copy them from an earlier, unrestricted
   Revision of dcraw.c, or (d) purchase a license from the author.

   The functions that process Foveon images have been RESTRICTED
   since Revision 1.237.  All other code remains free for all uses.

   *If you have not modified dcraw.c in any way, a link to my
   homepage qualifies as "full source code".

   $Revision: 1.444 $
   $Date: 2011/07/23 20:33:32 $
 */

#define DCRAW_VERSION "9.10"

#ifndef _GNU_SOURCE
#define _GNU_SOURCE
#endif
#define _USE_MATH_DEFINES
#include <ctype.h>
#include <errno.h>
#include <fcntl.h>
#include <float.h>
#include <limits.h>
#include <math.h>
#include <setjmp.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <sys/types.h>
#include <assert.h>
#define NO_JASPER
#define NO_LCMS

#ifdef NODEPS
#define NO_JASPER
#define NO_JPEG
#define NO_LCMS
#endif
#ifdef LOCALEDIR
#include <libintl.h>
#define _(String) gettext(String)
#else
#define _(String) (String)
#endif

#if defined(DJGPP) || defined(__MINGW32__)
#define fseeko fseek
#define ftello ftell
#else
#define fgetc getc_unlocked
#endif
#ifdef __CYGWIN__
#include <io.h>
#endif
#ifdef WIN32
#include <sys/utime.h>
#include <winsock2.h>
#pragma comment(lib, "ws2_32.lib")
#define snprintf _snprintf
#define strcasecmp stricmp
#define strncasecmp strnicmp
typedef __int64 INT64;
typedef unsigned __int64 UINT64;
#else
#include <unistd.h>
#include <utime.h>
#include <netinet/in.h>
typedef long long INT64;
typedef unsigned long long UINT64;
#endif

#ifdef LJPEG_DECODE
#error Please compile dcraw.c by itself.
#error Do not link it with ljpeg_decode.
#endif

#ifndef LONG_BIT
#define LONG_BIT (8 * sizeof (long))
#endif

#if !defined(uchar)
#define uchar unsigned char
#endif
#if !defined(ushort)
#define ushort unsigned short
#endif

/*
   All global variables are defined here, and all functions that
   access them are prefixed with "CLASS".  Note that a thread-safe
   C++ class cannot have non-const static local variables.
 */
//static FILE *ifp;
static short order;
static const char *ifname;
static char *meta_data;
static char cdesc[5], desc[512], make[64], model[64], model2[64], artist[64];
static float iso_speed, shutter, aperture, focal_len;
static time_t timestamp;
static unsigned shot_order, kodak_cbpp, filters, exif_cfa, unique_id;
static off_t    data_offset;
static off_t    thumb_offset, meta_offset, profile_offset;
static unsigned thumb_length, meta_length, profile_length;
static unsigned thumb_misc, fuji_layout;
static unsigned tiff_nifds, tiff_samples, tiff_bps, tiff_compress;
static unsigned black, cblack[8], maximum, mix_green, raw_color, zero_is_bad;
static unsigned zero_after_ff, is_raw, dng_version, is_foveon, data_error;
static unsigned tile_width, tile_length, gpsdata[32], load_flags;
static ushort raw_height, raw_width, height, width, top_margin, left_margin;
static ushort shrink, iheight, iwidth, fuji_width, thumb_width, thumb_height;
static int tiff_flip, colors;
static double pixel_aspect;
static ushort (*image)[4], white[8][8], curve[0x10000], cr2_slice[3], sraw_mul[4];
static int use_camera_wb=0;
static float cam_mul[4], pre_mul[4], cmatrix[3][4], rgb_cam[3][4];
static const double xyz_rgb[3][3] = {			/* XYZ from RGB */
  { 0.412453, 0.357580, 0.180423 },
  { 0.212671, 0.715160, 0.072169 },
  { 0.019334, 0.119193, 0.950227 } };
static const float d65_white[3] = { 0.950456, 1, 1.088754 };
static void (*load_raw)(ASTAR_FILE_RAM &fp);
static jmp_buf failure;

static struct tiff_ifd {
  int width, height, bps, comp, phint, offset, flip, samples, bytes;
} tiff_ifd[10];

#define CLASS

#define FORC(cnt) for (c=0; c < cnt; c++)
#define FORC3 FORC(3)
#define FORC4 FORC(4)
#define FORCC FORC(colors)

#define SQR(x) ((x)*(x))
#define SWAP(a,b) { a=a+b; b=a-b; a=a-b; }

/*
   In order to inline this calculation, I make the risky
   assumption that all filter patterns can be described
   by a repeating pattern of eight rows and two columns

   Do not use the FC or BAYER macros with the Leaf CatchLight,
   because its pattern is 16x16, not 2x8.

   Return values are either 0/1/2/3 = G/M/C/Y or 0/1/2/3 = R/G1/B/G2

	PowerShot 600	PowerShot A50	PowerShot Pro70	Pro90 & G1
	0xe1e4e1e4:	0x1b4e4b1e:	0x1e4b4e1b:	0xb4b4b4b4:

	  0 1 2 3 4 5	  0 1 2 3 4 5	  0 1 2 3 4 5	  0 1 2 3 4 5
	0 G M G M G M	0 C Y C Y C Y	0 Y C Y C Y C	0 G M G M G M
	1 C Y C Y C Y	1 M G M G M G	1 M G M G M G	1 Y C Y C Y C
	2 M G M G M G	2 Y C Y C Y C	2 C Y C Y C Y
	3 C Y C Y C Y	3 G M G M G M	3 G M G M G M
			4 C Y C Y C Y	4 Y C Y C Y C
	PowerShot A5	5 G M G M G M	5 G M G M G M
	0x1e4e1e4e:	6 Y C Y C Y C	6 C Y C Y C Y
			7 M G M G M G	7 M G M G M G
	  0 1 2 3 4 5
	0 C Y C Y C Y
	1 G M G M G M
	2 C Y C Y C Y
	3 M G M G M G

   All RGB cameras use one of these Bayer grids:

	0x16161616:	0x61616161:	0x49494949:	0x94949494:

	  0 1 2 3 4 5	  0 1 2 3 4 5	  0 1 2 3 4 5	  0 1 2 3 4 5
	0 B G B G B G	0 G R G R G R	0 G B G B G B	0 R G R G R G
	1 G R G R G R	1 B G B G B G	1 R G R G R G	1 G B G B G B
	2 B G B G B G	2 G R G R G R	2 G B G B G B	2 R G R G R G
	3 G R G R G R	3 B G B G B G	3 R G R G R G	3 G B G B G B
 */

#define FC(row,col) \
	(filters >> ((((row) << 1 & 14) + ((col) & 1)) << 1) & 3)

#define BAYER(row,col) \
	image[((row) >> shrink)*iwidth + ((col) >> shrink)][FC(row,col)]

#define BAYER2(row,col) \
	image[((row) >> shrink)*iwidth + ((col) >> shrink)][fc(row,col)]

#ifndef __GLIBC__
char *my_memmem (char *haystack, size_t haystacklen,
	      char *needle, size_t needlelen)
{
  char *c;
  for (c = haystack; c <= haystack + haystacklen - needlelen; c++)
    if (!memcmp (c, needle, needlelen))
      return c;
  return 0;
}
#define memmem my_memmem
#endif

void CLASS merror (void *ptr, const char *where)
{
  if (ptr) return;
  fprintf (stderr,_("%s: Out of memory in %s\n"), ifname, where);
  longjmp (failure, 1);
}

void CLASS derror( ASTAR_FILE_RAM &fp )
{
  if (!data_error) {
    fprintf (stderr, "%s: ", ifname);
    if ( fp.eof())
      fprintf (stderr,_("Unexpected end of file\n"));
    else
      fprintf (stderr,_("Corrupt data near 0x%llx\n"), (INT64) fp.tell());
  }
  data_error++;
}

ushort CLASS sget2 (uchar *s)
{
  if (order == 0x4949)		/* "II" means little-endian */
    return s[0] | s[1] << 8;
  else				/* "MM" means big-endian */
    return s[0] << 8 | s[1];
}

ushort CLASS get2( ASTAR_FILE_RAM &fp )
{
  if ( order == 0x4949 )
  {
    return fp.get_le_short();
  }
  return fp.get_be_short();
}

unsigned CLASS sget4 (uchar *s)
{
  if (order == 0x4949)
    return s[0] | s[1] << 8 | s[2] << 16 | s[3] << 24;
  else
    return s[0] << 24 | s[1] << 16 | s[2] << 8 | s[3];
}
#define sget4(s) sget4((uchar *)s)

unsigned CLASS get4( ASTAR_FILE_RAM &fp )
{
  uchar str[4] = { 0xff,0xff,0xff,0xff };
  //my_fread (str, 1, 4, ifp);
  fp.read( str, 4 );
  return sget4(str);
}

unsigned CLASS getint ( ASTAR_FILE_RAM &fp, int type)
{
  return type == 3 ? get2(fp) : get4(fp);
}

float CLASS int_to_float (int i)
{
  union { int i; float f; } u;
  u.i = i;
  return u.f;
}

double CLASS getreal ( ASTAR_FILE_RAM &fp, int type)
{
  union { char c[8]; double d; } u;
  int i, rev;

  switch (type) {
    case 3: return (unsigned short) get2( fp );
    case 4: return (unsigned int) get4( fp );
    case 5:  u.d = (unsigned int) get4( fp );
      return u.d / (unsigned int) get4( fp );
    case 8: return (signed short) get2( fp );
    case 9: return (signed int) get4( fp );
    case 10: u.d = (signed int) get4( fp );
      return u.d / (signed int) get4( fp );
    case 11: return int_to_float (get4( fp ));
    case 12:
      rev = 7 * ((order == 0x4949) == (ntohs(0x1234) == 0x1234));
      for (i=0; i < 8; i++)
	u.c[i ^ rev] = fp.getchr();
      return u.d;
    default: return fp.getchr();
  }
}

void CLASS read_shorts (ASTAR_FILE_RAM &fp, ushort *pixel, int count)
{
  fp.read (pixel, 2* count);
 	// < count) derror(fp);
  if ((order == 0x4949) == (ntohs(0x1234) == 0x1234))
    swab (pixel, pixel, count*2);
}

/*
   getbits(-1) initializes the buffer
   getbits(n) where 0 <= n <= 25 returns an n-bit integer
 */
unsigned CLASS getbithuff (ASTAR_FILE_RAM &fp, int nbits, ushort *huff)
{
  static unsigned bitbuf=0;
  static int vbits=0, reset=0;
  unsigned c;

  if (nbits == -1)
    return bitbuf = vbits = reset = 0;
  if (nbits == 0 || vbits < 0) return 0;
  while (!reset && vbits < nbits && (c = fp.getchr()) != EOF &&
    !(reset = zero_after_ff && c == 0xff && fp.getchr())) {
    bitbuf = (bitbuf << 8) + (uchar) c;
    vbits += 8;
  }
  c = bitbuf << (32-vbits) >> (32-nbits);
  if (huff) {
    vbits -= huff[c] >> 8;
    c = (uchar) huff[c];
  } else
    vbits -= nbits;
  if (vbits < 0) derror(fp);
  return c;
}

#define getbits(fp,n) getbithuff(fp,n,0)
#define gethuff(fp,h) getbithuff(fp,*h,h+1)

/*
   Construct a decode tree according the specification in *source.
   The first 16 bytes specify how many codes should be 1-bit, 2-bit
   3-bit, etc.  Bytes after that are the leaf values.

   For example, if the source is

    { 0,1,4,2,3,1,2,0,0,0,0,0,0,0,0,0,
      0x04,0x03,0x05,0x06,0x02,0x07,0x01,0x08,0x09,0x00,0x0a,0x0b,0xff  },

   then the code is

	00		0x04
	010		0x03
	011		0x05
	100		0x06
	101		0x02
	1100		0x07
	1101		0x01
	11100		0x08
	11101		0x09
	11110		0x00
	111110		0x0a
	1111110		0x0b
	1111111		0xff
 */
ushort * CLASS make_decoder_ref (const uchar **source)
{
  int max, len, h, i, j;
  const uchar *count;
  ushort *huff;

  count = (*source += 16) - 17;
  for (max=16; max && !count[max]; max--);
  huff = (ushort *) calloc (1 + (1 << max), sizeof *huff);
  merror (huff, "make_decoder()");
  huff[0] = max;
  for (h=len=1; len <= max; len++)
    for (i=0; i < count[len]; i++, ++*source)
      for (j=0; j < 1 << (max-len); j++)
	if (h <= 1 << max)
	  huff[h++] = len << 8 | **source;
  return huff;
}

ushort * CLASS make_decoder (const uchar *source)
{
  return make_decoder_ref (&source);
}

/*
   Not a full implementation of Lossless JPEG, just
   enough to decode Canon, Kodak and Adobe DNG images.
 */
struct jhead {
  int bits, high, wide, clrs, sraw, psv, restart, vpred[6];
  ushort *huff[6], *free[4], *row;
};

int CLASS ljpeg_start (ASTAR_FILE_RAM &fp, struct jhead *jh, int info_only)
{
  int c, tag, len;
  uchar data[0x10000];
  const uchar *dp;

  memset (jh, 0, sizeof *jh);
  jh->restart = INT_MAX;
  fp.read (data, 2);
  if (data[1] != 0xd8) return 0;
  do {
    fp.read (data, 2* 2);
    tag =  data[0] << 8 | data[1];
    len = (data[2] << 8 | data[3]) - 2;
    if (tag <= 0xff00) return 0;
    fp.read(data, len);
    switch (tag) {
      case 0xffc3:
	jh->sraw = ((data[7] >> 4) * (data[7] & 15) - 1) & 3;
      case 0xffc0:
	jh->bits = data[0];
	jh->high = data[1] << 8 | data[2];
	jh->wide = data[3] << 8 | data[4];
	jh->clrs = data[5] + jh->sraw;
	if (len == 9 && !dng_version) fp.getchr();
	break;
      case 0xffc4:
	if (info_only) break;
	for (dp = data; dp < data+len && (c = *dp++) < 4; )
	  jh->free[c] = jh->huff[c] = make_decoder_ref (&dp);
	break;
      case 0xffda:
	jh->psv = data[1+data[0]*2];
	jh->bits -= data[3+data[0]*2] & 15;
	break;
      case 0xffdd:
	jh->restart = data[0] << 8 | data[1];
    }
  } while (tag != 0xffda);
  if (info_only) return 1;
  FORC(5) if (!jh->huff[c+1]) jh->huff[c+1] = jh->huff[c];
  if (jh->sraw) {
    FORC(4)        jh->huff[2+c] = jh->huff[1];
    FORC(jh->sraw) jh->huff[1+c] = jh->huff[0];
  }
  jh->row = (ushort *) calloc (jh->wide*jh->clrs, 4);
  merror (jh->row, "ljpeg_start()");
  return zero_after_ff = 1;
}

void CLASS ljpeg_end (struct jhead *jh)
{
  int c;
  FORC4 if (jh->free[c]) free (jh->free[c]);
  free (jh->row);
}

int CLASS ljpeg_diff (ASTAR_FILE_RAM &fp, ushort *huff)
{
  int len, diff;

  len = gethuff(fp,huff);
  if (len == 16 && (!dng_version || dng_version >= 0x1010000))
    return -32768;
  diff = getbits(fp, len);
  if ((diff & (1 << (len-1))) == 0)
    diff -= (1 << len) - 1;
  return diff;
}

ushort * CLASS ljpeg_row (ASTAR_FILE_RAM &fp, int jrow, struct jhead *jh)
{
  int col, c, diff, pred, spred=0;
  ushort mark=0, *row[3];

  if (jrow * jh->wide % jh->restart == 0) {
    FORC(6) jh->vpred[c] = 1 << (jh->bits-1);
    if (jrow) {
      fp.seek (fp.tell()-2);
      do mark = (mark << 8) + (c = fp.getchr());
      while (c != EOF && mark >> 4 != 0xffd);
    }
    getbits(fp, -1);
  }
  FORC3 row[c] = jh->row + jh->wide*jh->clrs*((jrow+c) & 1);
  for (col=0; col < jh->wide; col++)
    FORC(jh->clrs) {
      diff = ljpeg_diff (fp, jh->huff[c]);
      if (jh->sraw && c <= jh->sraw && (col | c))
		    pred = spred;
      else if (col) pred = row[0][-jh->clrs];
      else	    pred = (jh->vpred[c] += diff) - diff;
      if (jrow && col) switch (jh->psv) {
	case 1:	break;
	case 2: pred = row[1][0];					break;
	case 3: pred = row[1][-jh->clrs];				break;
	case 4: pred = pred +   row[1][0] - row[1][-jh->clrs];		break;
	case 5: pred = pred + ((row[1][0] - row[1][-jh->clrs]) >> 1);	break;
	case 6: pred = row[1][0] + ((pred - row[1][-jh->clrs]) >> 1);	break;
	case 7: pred = (pred + row[1][0]) >> 1;				break;
	default: pred = 0;
      }
      if ((**row = pred + diff) >> jh->bits) derror(fp);
      if (c <= jh->sraw) spred = **row;
      row[0]++; row[1]++;
    }
  return row[2];
}

void CLASS lossless_jpeg_load_raw( ASTAR_FILE_RAM &fp )
{
  int jwide, jrow, jcol, val, jidx, c, i, j, row=0, col=0;
  struct jhead jh;
  int min=INT_MAX;
  ushort *rp;

  if (!ljpeg_start (fp, &jh, 0)) return;
  jwide = jh.wide * jh.clrs;

  for (jrow=0; jrow < jh.high; jrow++) {
    rp = ljpeg_row (fp, jrow, &jh);
    if (load_flags & 1)
      row = jrow & 1 ? height-1-jrow/2 : jrow/2;
    for (jcol=0; jcol < jwide; jcol++) {
      val = *rp++;
      if (jh.bits <= 12)
	val = curve[val & 0xfff];
      if (cr2_slice[0]) {
	jidx = jrow*jwide + jcol;
	i = jidx / (cr2_slice[1]*jh.high);
	if ((j = i >= cr2_slice[0]))
		 i  = cr2_slice[0];
	jidx -= i * (cr2_slice[1]*jh.high);
	row = jidx / cr2_slice[1+j];
	col = jidx % cr2_slice[1+j] + i*cr2_slice[1];
      }
      if (raw_width == 3984 && (col -= 2) < 0)
	col += (row--,raw_width);
      if ((unsigned) (row-top_margin) < height) {
	c = FC(row-top_margin,col-left_margin);
	if ((unsigned) (col-left_margin) < width) {
	  BAYER(row-top_margin,col-left_margin) = val;
	  if (min > val) min = val;
	} else if (col > 1 && (unsigned) (col-left_margin+2) > width+3)
	  cblack[c] += (cblack[4+c]++,val);
      }
      if (++col >= raw_width)
	col = (row++,0);
    }
  }
  ljpeg_end (&jh);
  FORC4 if (cblack[4+c]) cblack[c] /= cblack[4+c];
  if (!strcasecmp(make,"KODAK"))
    black = min;
}

void CLASS pseudoinverse (double (*in)[3], double (*out)[3], int size)
{
  double work[3][6], num;
  int i, j, k;

  for (i=0; i < 3; i++) {
    for (j=0; j < 6; j++)
      work[i][j] = j == i+3;
    for (j=0; j < 3; j++)
      for (k=0; k < size; k++)
	work[i][j] += in[k][i] * in[k][j];
  }
  for (i=0; i < 3; i++) {
    num = work[i][i];
    for (j=0; j < 6; j++)
      work[i][j] /= num;
    for (k=0; k < 3; k++) {
      if (k==i) continue;
      num = work[k][i];
      for (j=0; j < 6; j++)
	work[k][j] -= work[i][j] * num;
    }
  }
  for (i=0; i < size; i++)
    for (j=0; j < 3; j++)
      for (out[i][j]=k=0; k < 3; k++)
	out[i][j] += work[j][k+3] * in[i][k];
}

void CLASS cam_xyz_coeff (double cam_xyz[4][3])
{
  double cam_rgb[4][3], inverse[4][3], num;
  int i, j, k;

  for (i=0; i < colors; i++)		/* Multiply out XYZ colorspace */
    for (j=0; j < 3; j++)
      for (cam_rgb[i][j] = k=0; k < 3; k++)
	cam_rgb[i][j] += cam_xyz[i][k] * xyz_rgb[k][j];

  for (i=0; i < colors; i++) {		/* Normalize cam_rgb so that */
    for (num=j=0; j < 3; j++)		/* cam_rgb * (1,1,1) is (1,1,1,1) */
      num += cam_rgb[i][j];
    for (j=0; j < 3; j++)
      cam_rgb[i][j] /= num;
    pre_mul[i] = 1 / num;
  }
  pseudoinverse (cam_rgb, inverse, colors);
  for (raw_color = i=0; i < 3; i++)
    for (j=0; j < colors; j++)
      rgb_cam[i][j] = inverse[j][i];
}

void CLASS scale_colors()
{
  unsigned bottom, right, size, row, col, ur, uc, i, x, y, c, sum[8];
  int val, dark, sat;
  double dsum[8], dmin, dmax;
  float scale_mul[4], fr, fc;
  ushort *img=0, *pix;

  float user_mul[4]={0,0,0,0};
  int use_auto_wb = 0;
  unsigned greybox[4] = { 0, 0, UINT_MAX, UINT_MAX };
#define MIN(a,b) ((a) < (b) ? (a) : (b))
#define MAX(a,b) ((a) > (b) ? (a) : (b))
#define LIM(x,min,max) MAX(min,MIN(x,max))

  int highlight=0;
  int verbose=0;
#define CLIP(x) LIM(x,0,65535)
  double aber[4]={1,1,1,1};

  FORC4 cblack[c] += black;
  if (user_mul[0])
    memcpy (pre_mul, user_mul, sizeof pre_mul);
  if (use_auto_wb || (use_camera_wb && cam_mul[0] == -1)) {
    memset (dsum, 0, sizeof dsum);
    bottom = MIN (greybox[1]+greybox[3], height);
    right  = MIN (greybox[0]+greybox[2], width);
    for (row=greybox[1]; row < bottom; row += 8)
      for (col=greybox[0]; col < right; col += 8) {
        memset (sum, 0, sizeof sum);
        for (y=row; y < row+8 && y < bottom; y++)
          for (x=col; x < col+8 && x < right; x++)
            FORC4 {
              if (filters) {
                c = FC(y,x);
                val = BAYER(y,x);
              } else
                val = image[y*width+x][c];
              if (val > maximum-25) goto skip_block;
              if ((val -= cblack[c]) < 0) val = 0;
              sum[c] += val;
              sum[c+4]++;
              if (filters) break;
            }
        FORC(8) dsum[c] += sum[c];
skip_block: ;
      }
    FORC4 if (dsum[c]) pre_mul[c] = dsum[c+4] / dsum[c];
  }
  if (use_camera_wb && cam_mul[0] != -1) {
    memset (sum, 0, sizeof sum);
    for (row=0; row < 8; row++)
      for (col=0; col < 8; col++) {
        c = FC(row,col);
        if ((val = white[row][col] - cblack[c]) > 0)
          sum[c] += val;
        sum[c+4]++;
      }
    if (sum[0] && sum[1] && sum[2] && sum[3])
      FORC4 pre_mul[c] = (float) sum[c+4] / sum[c];
    else if (cam_mul[0] && cam_mul[2])
      memcpy (pre_mul, cam_mul, sizeof pre_mul);
    else
      fprintf (stderr,_("%s: Cannot use camera white balance.\n"), ifname);
  }
  if (pre_mul[3] == 0) pre_mul[3] = colors < 4 ? pre_mul[1] : 1;
  dark = black;
  sat = maximum;
  maximum -= black;
  for (dmin=DBL_MAX, dmax=c=0; c < 4; c++) {
    if (dmin > pre_mul[c])
        dmin = pre_mul[c];
    if (dmax < pre_mul[c])
        dmax = pre_mul[c];
  }
  if (!highlight) dmax = dmin;
  FORC4 scale_mul[c] = (pre_mul[c] /= dmax) * 65535.0 / maximum;
  if (verbose) {
    fprintf (stderr,
      _("Scaling with darkness %d, saturation %d, and\nmultipliers"), dark, sat);
    FORC4 fprintf (stderr, " %f", pre_mul[c]);
    fputc ('\n', stderr);
  }
  size = iheight*iwidth;
  for (i=0; i < size*4; i++) {
    val = image[0][i];
    if (!val) continue;
    val -= cblack[i & 3];
    val *= scale_mul[i & 3];
    image[0][i] = CLIP(val);
  }
  if ((aber[0] != 1 || aber[2] != 1) && colors == 3) {
    if (verbose)
      fprintf (stderr,_("Correcting chromatic aberration...\n"));
    for (c=0; c < 4; c+=2) {
      if (aber[c] == 1) continue;
      img = (ushort *) malloc (size * sizeof *img);
      merror (img, "scale_colors()");
      for (i=0; i < size; i++)
        img[i] = image[i][c];
      for (row=0; row < iheight; row++) {
        ur = fr = (row - iheight*0.5) * aber[c] + iheight*0.5;
        if (ur > iheight-2) continue;
        fr -= ur;
        for (col=0; col < iwidth; col++) {
          uc = fc = (col - iwidth*0.5) * aber[c] + iwidth*0.5;
          if (uc > iwidth-2) continue;
          fc -= uc;
          pix = img + ur*iwidth + uc;
          image[row*iwidth+col][c] =
            (pix[     0]*(1-fc) + pix[       1]*fc) * (1-fr) +
            (pix[iwidth]*(1-fc) + pix[iwidth+1]*fc) * fr;
        }
      }
      free(img);
    }
  }
}


void CLASS tiff_get (
	ASTAR_FILE_RAM &fp,
	unsigned base,
	unsigned *tag, unsigned *type, unsigned *len, unsigned *save)
{
  *tag  = get2(fp);
  *type = get2(fp);
  *len  = get4(fp);
  *save = fp.tell() + 4;
  if (*len * ("11124811248488"[*type < 14 ? *type:0]-'0') > 4)
    fp.seek (get4(fp)+base);
}

void CLASS parse_thumb_note (ASTAR_FILE_RAM &fp, int base, unsigned toff, unsigned tlen)
{
  unsigned entries, tag, type, len, save;

  entries = get2(fp);
  while (entries--) {
    tiff_get (fp, base, &tag, &type, &len, &save);
    if (tag == toff) thumb_offset = get4(fp)+base;
    if (tag == tlen) thumb_length = get4(fp);
    fp.seek (save);
  }
}

int CLASS parse_tiff_ifd (ASTAR_FILE_RAM &fp, int base);

void CLASS parse_makernote (ASTAR_FILE_RAM &fp, int base, int uptag)
{
  static const uchar xlat[2][256] = {
  { 0xc1,0xbf,0x6d,0x0d,0x59,0xc5,0x13,0x9d,0x83,0x61,0x6b,0x4f,0xc7,0x7f,0x3d,0x3d,
    0x53,0x59,0xe3,0xc7,0xe9,0x2f,0x95,0xa7,0x95,0x1f,0xdf,0x7f,0x2b,0x29,0xc7,0x0d,
    0xdf,0x07,0xef,0x71,0x89,0x3d,0x13,0x3d,0x3b,0x13,0xfb,0x0d,0x89,0xc1,0x65,0x1f,
    0xb3,0x0d,0x6b,0x29,0xe3,0xfb,0xef,0xa3,0x6b,0x47,0x7f,0x95,0x35,0xa7,0x47,0x4f,
    0xc7,0xf1,0x59,0x95,0x35,0x11,0x29,0x61,0xf1,0x3d,0xb3,0x2b,0x0d,0x43,0x89,0xc1,
    0x9d,0x9d,0x89,0x65,0xf1,0xe9,0xdf,0xbf,0x3d,0x7f,0x53,0x97,0xe5,0xe9,0x95,0x17,
    0x1d,0x3d,0x8b,0xfb,0xc7,0xe3,0x67,0xa7,0x07,0xf1,0x71,0xa7,0x53,0xb5,0x29,0x89,
    0xe5,0x2b,0xa7,0x17,0x29,0xe9,0x4f,0xc5,0x65,0x6d,0x6b,0xef,0x0d,0x89,0x49,0x2f,
    0xb3,0x43,0x53,0x65,0x1d,0x49,0xa3,0x13,0x89,0x59,0xef,0x6b,0xef,0x65,0x1d,0x0b,
    0x59,0x13,0xe3,0x4f,0x9d,0xb3,0x29,0x43,0x2b,0x07,0x1d,0x95,0x59,0x59,0x47,0xfb,
    0xe5,0xe9,0x61,0x47,0x2f,0x35,0x7f,0x17,0x7f,0xef,0x7f,0x95,0x95,0x71,0xd3,0xa3,
    0x0b,0x71,0xa3,0xad,0x0b,0x3b,0xb5,0xfb,0xa3,0xbf,0x4f,0x83,0x1d,0xad,0xe9,0x2f,
    0x71,0x65,0xa3,0xe5,0x07,0x35,0x3d,0x0d,0xb5,0xe9,0xe5,0x47,0x3b,0x9d,0xef,0x35,
    0xa3,0xbf,0xb3,0xdf,0x53,0xd3,0x97,0x53,0x49,0x71,0x07,0x35,0x61,0x71,0x2f,0x43,
    0x2f,0x11,0xdf,0x17,0x97,0xfb,0x95,0x3b,0x7f,0x6b,0xd3,0x25,0xbf,0xad,0xc7,0xc5,
    0xc5,0xb5,0x8b,0xef,0x2f,0xd3,0x07,0x6b,0x25,0x49,0x95,0x25,0x49,0x6d,0x71,0xc7 },
  { 0xa7,0xbc,0xc9,0xad,0x91,0xdf,0x85,0xe5,0xd4,0x78,0xd5,0x17,0x46,0x7c,0x29,0x4c,
    0x4d,0x03,0xe9,0x25,0x68,0x11,0x86,0xb3,0xbd,0xf7,0x6f,0x61,0x22,0xa2,0x26,0x34,
    0x2a,0xbe,0x1e,0x46,0x14,0x68,0x9d,0x44,0x18,0xc2,0x40,0xf4,0x7e,0x5f,0x1b,0xad,
    0x0b,0x94,0xb6,0x67,0xb4,0x0b,0xe1,0xea,0x95,0x9c,0x66,0xdc,0xe7,0x5d,0x6c,0x05,
    0xda,0xd5,0xdf,0x7a,0xef,0xf6,0xdb,0x1f,0x82,0x4c,0xc0,0x68,0x47,0xa1,0xbd,0xee,
    0x39,0x50,0x56,0x4a,0xdd,0xdf,0xa5,0xf8,0xc6,0xda,0xca,0x90,0xca,0x01,0x42,0x9d,
    0x8b,0x0c,0x73,0x43,0x75,0x05,0x94,0xde,0x24,0xb3,0x80,0x34,0xe5,0x2c,0xdc,0x9b,
    0x3f,0xca,0x33,0x45,0xd0,0xdb,0x5f,0xf5,0x52,0xc3,0x21,0xda,0xe2,0x22,0x72,0x6b,
    0x3e,0xd0,0x5b,0xa8,0x87,0x8c,0x06,0x5d,0x0f,0xdd,0x09,0x19,0x93,0xd0,0xb9,0xfc,
    0x8b,0x0f,0x84,0x60,0x33,0x1c,0x9b,0x45,0xf1,0xf0,0xa3,0x94,0x3a,0x12,0x77,0x33,
    0x4d,0x44,0x78,0x28,0x3c,0x9e,0xfd,0x65,0x57,0x16,0x94,0x6b,0xfb,0x59,0xd0,0xc8,
    0x22,0x36,0xdb,0xd2,0x63,0x98,0x43,0xa1,0x04,0x87,0x86,0xf7,0xa6,0x26,0xbb,0xd6,
    0x59,0x4d,0xbf,0x6a,0x2e,0xaa,0x2b,0xef,0xe6,0x78,0xb6,0x4e,0xe0,0x2f,0xdc,0x7c,
    0xbe,0x57,0x19,0x32,0x7e,0x2a,0xd0,0xb8,0xba,0x29,0x00,0x3c,0x52,0x7d,0xa8,0x49,
    0x3b,0x2d,0xeb,0x25,0x49,0xfa,0xa3,0xaa,0x39,0xa7,0xc5,0xa7,0x50,0x11,0x36,0xfb,
    0xc6,0x67,0x4a,0xf5,0xa5,0x12,0x65,0x7e,0xb0,0xdf,0xaf,0x4e,0xb3,0x61,0x7f,0x2f } };
  unsigned offset=0, entries, tag, type, len, save, c;
  unsigned ver97=0, serial=0, i, wbi=0, wb[4]={0,0,0,0};
  uchar buf97[324], ci, cj, ck;
  short morder, sorder=order;
  char buf[10];
/*
   The MakerNote might have its own TIFF header (possibly with
   its own byte-order!), or it might just be a table.
 */
  fp.read(buf, 10);
  if (!strncmp (buf,"KDK" ,3) ||	/* these aren't TIFF tables */
      !strncmp (buf,"VER" ,3) ||
      !strncmp (buf,"IIII",4) ||
      !strncmp (buf,"MMMM",4)) return;
  if (!strncmp (buf,"KC"  ,2) ||	/* Konica KD-400Z, KD-510Z */
      !strncmp (buf,"MLY" ,3)) {	/* Minolta DiMAGE G series */
    order = 0x4d4d;
    while ((i=fp.tell()) < data_offset && i < 16384) {
      wb[0] = wb[2];  wb[2] = wb[1];  wb[1] = wb[3];
      wb[3] = get2(fp);
      if (wb[1] == 256 && wb[3] == 256 &&
	  wb[0] > 256 && wb[0] < 640 && wb[2] > 256 && wb[2] < 640)
	FORC4 cam_mul[c] = wb[c];
    }
    goto quit;
  }
  if (!strcmp (buf,"Nikon")) {
    base = fp.tell();
    order = get2(fp);
    if (get2(fp) != 42) goto quit;
    offset = get4(fp);
    fp.seek (fp.tell() + offset-8);
  } else if (!strcmp (buf,"OLYMPUS")) {
    base = fp.tell()-10;
    fp.seek( fp.tell() -2);
    order = get2(fp);  get2(fp);
  } else if (!strncmp (buf,"SONY",4) ||
	     !strcmp  (buf,"Panasonic")) {
    goto nf;
  } else if (!strncmp (buf,"FUJIFILM",8)) {
    base = fp.tell()-10;
nf: order = 0x4949;
    fp.seek (fp.tell() + 2);
  } else if (!strcmp (buf,"OLYMP") ||
	     !strcmp (buf,"LEICA") ||
	     !strcmp (buf,"Ricoh") ||
	     !strcmp (buf,"EPSON"))
    fp.seek (fp.tell() -2);
  else if (!strcmp (buf,"AOC") ||
	   !strcmp (buf,"QVC"))
    fp.seek (fp.tell() -4);
  else {
    fp.seek (fp.tell() -10);
    if (!strncmp(make,"SAMSUNG",7))
      base = fp.tell();
  }
  entries = get2(fp);
  if (entries > 1000) return;
  morder = order;
  while (entries--) {
    order = morder;
    tiff_get (fp, base, &tag, &type, &len, &save);
    tag |= uptag << 16;
    if (tag == 2 && strstr(make,"NIKON") && !iso_speed)
      iso_speed = (get2(fp),get2(fp));
    if (tag == 4 && len > 26 && len < 35) {
      if ((i=(get4(fp),get2(fp))) != 0x7fff && !iso_speed)
	iso_speed = 50 * pow (2, i/32.0 - 4);
      if ((i=(get2(fp),get2(fp))) != 0x7fff && !aperture)
	aperture = pow (2, i/64.0);
      if ((i=get2(fp)) != 0xffff && !shutter)
	shutter = pow (2, (short) i/-32.0);
      wbi = (get2(fp),get2(fp));
      shot_order = (get2(fp),get2(fp));
    }
    if (tag == 7 && type == 2 && len > 20)
      fp.gets (model2, 64);
    if (tag == 8 && type == 4)
      shot_order = get4(fp);
    if (tag == 9 && !strcmp(make,"Canon"))
      fp.read (artist, 64);
    if (tag == 0xc && len == 4) {
      cam_mul[0] = getreal(fp,type);
      cam_mul[2] = getreal(fp,type);
    }
    if (tag == 0x10 && type == 4)
      unique_id = get4(fp);
    if (tag == 0x11 && is_raw && !strncmp(make,"NIKON",5)) {
      fp.seek (get4(fp)+base);
      parse_tiff_ifd (fp, base);
    }
    if (tag == 0x14 && type == 7) {
      if (len == 2560) {
	fp.seek (fp.tell() + 1248);
	goto get2_256;
      }
      fp.read (buf, 10);
      if (!strncmp(buf,"NRW ",4)) {
	fp.seek (fp.tell() + (strcmp(buf+4,"0100") ? 46:1546));
	cam_mul[0] = get4(fp) << 2;
	cam_mul[1] = get4(fp) + get4(fp);
	cam_mul[2] = get4(fp) << 2;
      }
    }
    if (tag == 0x15 && type == 2 && is_raw)
      fp.read (model, 64);
    if (strstr(make,"PENTAX")) {
      if (tag == 0x1b) tag = 0x1018;
      if (tag == 0x1c) tag = 0x1017;
    }
    if (tag == 0x1d)
      while ((c = fp.getchr()) && c != EOF)
	serial = serial*10 + (isdigit(c) ? c - '0' : c % 10);
    if (tag == 0x81 && type == 4) {
      data_offset = get4(fp);
      fp.seek (data_offset + 41);
      raw_height = get2(fp) * 2;
      raw_width  = get2(fp);
      filters = 0x61616161;
    }
    if (tag == 0x29 && type == 1) {
      c = wbi < 18 ? "012347800000005896"[wbi]-'0' : 0;
      fp.seek (fp.tell() + 8 + c*32);
      FORC4 cam_mul[c ^ (c >> 1) ^ 1] = get4(fp);
    }
    if ((tag == 0x81  && type == 7) ||
	(tag == 0x100 && type == 7) ||
	(tag == 0x280 && type == 1)) {
      thumb_offset = fp.tell();
      thumb_length = len;
    }
    if (tag == 0x88 && type == 4 && (thumb_offset = get4(fp)))
      thumb_offset += base;
    if (tag == 0x89 && type == 4)
      thumb_length = get4(fp);
    if (tag == 0x8c || tag == 0x96)
      meta_offset = fp.tell();
    if (tag == 0x97) {
      for (i=0; i < 4; i++)
	ver97 = ver97 * 10 + fp.getchr()-'0';
      switch (ver97) {
	case 100:
	  fp.seek (fp.tell() + 68);
	  FORC4 cam_mul[(c >> 1) | ((c & 1) << 1)] = get2(fp);
	  break;
	case 102:
	  fp.seek (fp.tell() + 6);
	  goto get2_rggb;
	case 103:
	  fp.seek (fp.tell() + 16);
	  FORC4 cam_mul[c] = get2(fp);
      }
      if (ver97 >= 200) {
	if (ver97 != 205) fp.seek (fp.tell() + 280);
	fp.read (buf97, 324);
      }
    }
    if (tag == 0xa1 && type == 7) {
      order = 0x4949;
      fp.seek (fp.tell() + 140);
      FORC3 cam_mul[c] = get4(fp);
    }
    if (tag == 0xa4 && type == 3) {
      fp.seek (fp.tell() + wbi*48);
      FORC3 cam_mul[c] = get2(fp);
    }
    if (tag == 0xa7 && (unsigned) (ver97-200) < 17) {
      ci = xlat[0][serial & 0xff];
      cj = xlat[1][fp.getchr()^fp.getchr()^fp.getchr()^fp.getchr()];
      ck = 0x60;
      for (i=0; i < 324; i++)
	buf97[i] ^= (cj += ci * ck++);
      i = "66666>666;6A;:;55"[ver97-200] - '0';
      FORC4 cam_mul[c ^ (c >> 1) ^ (i & 1)] =
	sget2 (buf97 + (i & -2) + c*2);
    }
    if (tag == 0x200 && len == 3)
      shot_order = (get4(fp),get4(fp));
    if (tag == 0x200 && len == 4)
      FORC4 cblack[c ^ c >> 1] = get2(fp);
    if (tag == 0x201 && len == 4)
      goto get2_rggb;
    if (tag == 0x220 && type == 7)
      meta_offset = fp.tell();
    if (tag == 0x401 && type == 4 && len == 4)
      FORC4 cblack[c ^ c >> 1] = get4(fp);
    if (tag == 0xe80 && len == 256 && type == 7) {
      fp.seek (fp.tell() + 48);
      cam_mul[0] = get2(fp) * 508 * 1.078 / 0x10000;
      cam_mul[2] = get2(fp) * 382 * 1.173 / 0x10000;
    }
    if (tag == 0xf00 && type == 7) {
      if (len == 614)
	fp.seek (fp.tell() + 176);
      else if (len == 734 || len == 1502)
	fp.seek (fp.tell() + 148);
      else goto next;
      goto get2_256;
    }
    if ((tag == 0x1011 && len == 9) || tag == 0x20400200)
      for (i=0; i < 3; i++)
	FORC3 cmatrix[i][c] = ((short) get2(fp)) / 256.0;
    if ((tag == 0x1012 || tag == 0x20400600) && len == 4)
      FORC4 cblack[c ^ c >> 1] = get2(fp);
    if (tag == 0x1017 || tag == 0x20400100)
      cam_mul[0] = get2(fp) / 256.0;
    if (tag == 0x1018 || tag == 0x20400100)
      cam_mul[2] = get2(fp) / 256.0;
    if (tag == 0x2011 && len == 2) {
get2_256:
      order = 0x4d4d;
      cam_mul[0] = get2(fp) / 256.0;
      cam_mul[2] = get2(fp) / 256.0;
    }
    if ((tag | 0x70) == 0x2070 && type == 4)
      fp.seek (get4(fp)+base);
    if (tag == 0x2020)
      parse_thumb_note (fp, base, 257, 258);
    if (tag == 0x2040)
      parse_makernote (fp, base, 0x2040);
    if (tag == 0xb028) {
      fp.seek (get4(fp)+base);
      parse_thumb_note (fp, base, 136, 137);
    }
    if (tag == 0x4001 && len > 500) {
      i = len == 582 ? 50 : len == 653 ? 68 : len == 5120 ? 142 : 126;
      fp.seek (i + fp.tell());
get2_rggb:
      FORC4 cam_mul[c ^ (c >> 1)] = get2(fp);
      fp.seek (22 + fp.tell());
      FORC4 sraw_mul[c ^ (c >> 1)] = get2(fp);
    }
    if (tag == 0xa021)
      FORC4 cam_mul[c ^ (c >> 1)] = get4(fp);
    if (tag == 0xa028)
      FORC4 cam_mul[c ^ (c >> 1)] -= get4(fp);
next:
    fp.seek (save);
  }
quit:
  order = sorder;
}

/*
   Since the TIFF DateTime string has no timezone information,
   assume that the camera's clock was set to Universal Time.
 */
void CLASS get_timestamp (ASTAR_FILE_RAM &fp, int reversed)
{
  struct tm t;
  char str[20];
  int i;

  str[19] = 0;
  if (reversed)
    for (i=19; i--; ) str[i] = fp.getchr();
  else
    fp.read (str, 19);
  memset (&t, 0, sizeof t);
  if (sscanf (str, "%d:%d:%d %d:%d:%d", &t.tm_year, &t.tm_mon,
	&t.tm_mday, &t.tm_hour, &t.tm_min, &t.tm_sec) != 6)
    return;
  t.tm_year -= 1900;
  t.tm_mon -= 1;
  t.tm_isdst = -1;
  if (mktime(&t) > 0)
    timestamp = mktime(&t);
}

void CLASS parse_exif ( ASTAR_FILE_RAM &fp, int base)
{
  unsigned kodak, entries, tag, type, len, save, c;
  double expo;

  kodak = !strncmp(make,"EASTMAN",7) && tiff_nifds < 3;
  entries = get2(fp);
  while (entries--) {
    tiff_get (fp, base, &tag, &type, &len, &save);
    switch (tag) {
      case 33434:  shutter = getreal(fp,type);		break;
      case 33437:  aperture = getreal(fp,type);		break;
      case 34855:  iso_speed = get2(fp);			break;
      case 36867:
      case 36868:  get_timestamp(fp, 0);			break;
      case 37377:  if ((expo = -getreal(fp,type)) < 128)
		     shutter = pow (2, expo);		break;
      case 37378:  aperture = pow (2, getreal(fp,type)/2);	break;
      case 37386:  focal_len = getreal(fp,type);		break;
      case 37500:  parse_makernote (fp, base, 0);		break;
      case 40962:  if (kodak) raw_width  = get4(fp);	break;
      case 40963:  if (kodak) raw_height = get4(fp);	break;
      case 41730:
	if (get4(fp) == 0x20002)
	  for (exif_cfa=c=0; c < 8; c+=2)
	    exif_cfa |= fp.getchr() * 0x01010101 << c;
    }
    fp.seek (save);
  }
}

int CLASS parse_tiff (ASTAR_FILE_RAM &fp, int base);

int CLASS parse_tiff_ifd (ASTAR_FILE_RAM &fp, int base)
{
  unsigned entries, tag, type, len, plen=16, save;
  int ifd, use_cm=0, cfa, i, j, c, ima_len=0;
  int blrr=1, blrc=1, dblack[] = { 0,0,0,0 };
  char software[64], *cbuf, *cp;
  uchar cfa_pat[16], cfa_pc[] = { 0,1,2,3 }, tab[256];
  double cc[4][4], cm[4][3], cam_xyz[4][3], num;
  double ab[]={ 1,1,1,1 }, asn[] = { 0,0,0,0 }, xyz[] = { 1,1,1 };
  unsigned sony_curve[] = { 0,0,0,0,0,4095 };
  unsigned sony_offset, sony_length, sony_key;
  struct jhead jh;
  FILE *sfp;

  if (tiff_nifds >= sizeof tiff_ifd / sizeof tiff_ifd[0])
    return 1;
  ifd = tiff_nifds++;
  for (j=0; j < 4; j++)
    for (i=0; i < 4; i++)
      cc[j][i] = i == j;
  entries = get2(fp);
  if (entries > 512) return 1;
  while (entries--) {
    tiff_get (fp, base, &tag, &type, &len, &save);
    switch (tag) {
      case 5:   width  = get2( fp );  break;
      case 6:   height = get2( fp );  break;
      case 7:   width += get2( fp );  break;
      case 9:  filters = get2( fp );  break;
      case 17: case 18:
	if (type == 3 && len == 1)
	  cam_mul[(tag-17)*2] = get2( fp ) / 256.0;
	break;
      case 23:
	if (type == 3) iso_speed = get2( fp );
	break;
      case 36: case 37: case 38:
	cam_mul[tag-0x24] = get2( fp );
	break;
      case 39:
	if (len < 50 || cam_mul[0]) break;
	fp.seek (12 + fp.tell());
	FORC3 cam_mul[c] = get2( fp );
	break;
      case 46:
	if (type != 7 || fp.getchr() != 0xff || fp.getchr() != 0xd8) break;
	thumb_offset = fp.tell() - 2;
	thumb_length = len;
	break;
      case 61440:			/* Fuji HS10 table */
	parse_tiff_ifd (fp, base);
	break;
      case 2: case 256: case 61441:	/* ImageWidth */
	tiff_ifd[ifd].width = getint(fp,type);
	break;
      case 3: case 257: case 61442:	/* ImageHeight */
	tiff_ifd[ifd].height = getint(fp,type);
	break;
      case 258:				/* BitsPerSample */
      case 61443:
	tiff_ifd[ifd].samples = len & 7;
	tiff_ifd[ifd].bps = getint(fp,type);
	break;
      case 259:				/* Compression */
	tiff_ifd[ifd].comp = get2(fp);
	break;
      case 262:				/* PhotometricInterpretation */
	tiff_ifd[ifd].phint = get2(fp);
	break;
      case 270:				/* ImageDescription */
	fp.read (desc, 512);
	break;
      case 271:				/* Make */
	fp.gets (make, 64);
	break;
      case 272:				/* Model */
	fp.gets (model, 64);
	break;
      case 273:				/* StripOffset */
      case 513:				/* JpegIFOffset */
      case 61447:
	tiff_ifd[ifd].offset = get4(fp)+base;
	if (!tiff_ifd[ifd].bps && tiff_ifd[ifd].offset > 0) {
	  fp.seek (tiff_ifd[ifd].offset);
	  if (ljpeg_start (fp, &jh, 1)) {
	    tiff_ifd[ifd].comp    = 6;
	    tiff_ifd[ifd].width   = jh.wide;
	    tiff_ifd[ifd].height  = jh.high;
	    tiff_ifd[ifd].bps     = jh.bits;
	    tiff_ifd[ifd].samples = jh.clrs;
	    if (!(jh.sraw || (jh.clrs & 1)))
	      tiff_ifd[ifd].width *= jh.clrs;
	    i = order;
	    parse_tiff (fp, tiff_ifd[ifd].offset + 12);
	    order = i;
	  }
	}
	break;
      case 274:				/* Orientation */
	tiff_ifd[ifd].flip = "50132467"[get2(fp) & 7]-'0';
	break;
      case 277:				/* SamplesPerPixel */
	tiff_ifd[ifd].samples = getint(fp,type) & 7;
	break;
      case 279:				/* StripByteCounts */
      case 514:
      case 61448:
	tiff_ifd[ifd].bytes = get4(fp);
	break;
      case 61454:
	FORC3 cam_mul[(4-c) % 3] = getint(fp,type);
	break;
      case 305:  case 11:		/* Software */
	fp.gets (software, 64);
	if (!strncmp(software,"Adobe",5) ||
	    !strncmp(software,"dcraw",5) ||
	    !strncmp(software,"UFRaw",5) ||
	    !strncmp(software,"Bibble",6) ||
	    !strncmp(software,"Nikon Scan",10) ||
	    !strcmp (software,"Digital Photo Professional"))
	  is_raw = 0;
	break;
      case 306:				/* DateTime */
	get_timestamp(fp, 0);
	break;
      case 315:				/* Artist */
	fp.read(artist, 64);
	break;
      case 322:				/* TileWidth */
	tile_width = getint(fp,type);
	break;
      case 323:				/* TileLength */
	tile_length = getint(fp,type);
	break;
      case 330:				/* SubIFDs */
	while (len--) {
	  i = fp.tell();
	  fp.seek (get4(fp)+base);
	  if (parse_tiff_ifd (fp, base)) break;
	  fp.seek (i+4);
	}
	break;
      case 400:
	strcpy (make, "Sarnoff");
	maximum = 0xfff;
	break;
      case 28688:
	FORC4 sony_curve[c+1] = get2(fp) >> 2 & 0xfff;
	for (i=0; i < 5; i++)
	  for (j = sony_curve[i]+1; j <= sony_curve[i+1]; j++)
	    curve[j] = curve[j-1] + (1 << i);
	break;
      case 29184: sony_offset = get4(fp);  break;
      case 29185: sony_length = get4(fp);  break;
      case 29217: sony_key    = get4(fp);  break;
      case 29443:
	FORC4 cam_mul[c ^ (c < 2)] = get2(fp);
	break;
      case 29459:
	FORC4 cam_mul[c] = get2(fp);
	i = (cam_mul[1] == 1024 && cam_mul[2] == 1024) << 1;
	SWAP (cam_mul[i],cam_mul[i+1])
	break;
      case 33405:			/* Model2 */
	fp.gets (model2, 64);
	break;
      case 33422:			/* CFAPattern */
      case 64777:			/* Kodak P-series */
	if ((plen=len) > 16) plen = 16;
	fp.read(cfa_pat, plen);
	for (colors=cfa=i=0; i < plen; i++) {
	  colors += !(cfa & (1 << cfa_pat[i]));
	  cfa |= 1 << cfa_pat[i];
	}
	if (cfa == 070) memcpy (cfa_pc,"\003\004\005",3);	/* CMY */
	if (cfa == 072) memcpy (cfa_pc,"\005\003\004\001",4);	/* GMCY */
	goto guess_cfa_pc;
      case 33424:
      case 33434:			/* ExposureTime */
	shutter = getreal(fp,type);
	break;
      case 33437:			/* FNumber */
	aperture = getreal(fp,type);
	break;
      case 34306:			/* Leaf white balance */
	FORC4 cam_mul[c ^ 1] = 4096.0 / get2(fp);
	break;
      case 34307:			/* Leaf CatchLight color matrix */
	fp.read (software, 7);
	if (strncmp(software,"MATRIX",6)) break;
	colors = 4;
	for (raw_color = i=0; i < 3; i++) {
	  //FORC4 fscanf (ifp, "%f", &rgb_cam[i][c^1]);
	  FORC4 rgb_cam[i][c^1]=fp.getfloat();
	  if (!use_camera_wb) continue;
	  num = 0;
	  FORC4 num += rgb_cam[i][c];
	  FORC4 rgb_cam[i][c] /= num;
	}
	break;
      case 34303:
	strcpy (make, "Leaf");
	break;
      case 34665:			/* EXIF tag */
	fp.seek (get4(fp)+base);
	parse_exif (fp, base);
	break;
      case 34675:			/* InterColorProfile */
      case 50831:			/* AsShotICCProfile */
	profile_offset = fp.tell();
	profile_length = len;
	break;
      case 37122:			/* CompressedBitsPerPixel */
	kodak_cbpp = get4(fp);
	break;
      case 37386:			/* FocalLength */
	focal_len = getreal(fp,type);
	break;
      case 37393:			/* ImageNumber */
	shot_order = getint(fp,type);
	break;
      case 37400:			/* old Kodak KDC tag */
	for (raw_color = i=0; i < 3; i++) {
	  getreal(fp,type);
	  FORC3 rgb_cam[i][c] = getreal(fp,type);
	}
	break;
      case 46275:			/* Imacon tags */
	strcpy (make, "Imacon");
	data_offset = fp.tell();
	ima_len = len;
	break;
      case 46279:
        assert( 0 );
	break;
      case 50454:			/* Sinar tag */
      case 50455:
	if (!(cbuf = (char *) malloc(len))) break;
	fp.read(cbuf, len);
	for (cp = cbuf-1; cp && cp < cbuf+len; cp = strchr(cp,'\n'))
	  if (!strncmp (++cp,"Neutral ",8))
	    sscanf (cp+8, "%f %f %f", cam_mul, cam_mul+1, cam_mul+2);
	free (cbuf);
	break;
      case 50458:
	if (!make[0]) strcpy (make, "Hasselblad");
	break;
      case 50459:			/* Hasselblad tag */
	i = order;
	j = fp.tell();
	c = tiff_nifds;
	order = get2(fp);
	fp.seek(j+(get2(fp),get4(fp)));
	parse_tiff_ifd (fp, j);
	maximum = 0xffff;
	tiff_nifds = c;
	order = i;
	break;
      case 50706:			/* DNGVersion */
	FORC4 dng_version = (dng_version << 8) + fp.getchr();
	if (!make[0]) strcpy (make, "DNG");
	is_raw = 1;
	break;
      case 50710:			/* CFAPlaneColor */
	if (len > 4) len = 4;
	colors = len;
	fp.read (cfa_pc, colors);
guess_cfa_pc:
	FORCC tab[cfa_pc[c]] = c;
	cdesc[c] = 0;
	for (i=16; i--; )
	  filters = filters << 2 | tab[cfa_pat[i % plen]];
	break;
      case 50711:			/* CFALayout */
	if (get2( fp ) == 2) {
	  fuji_width = 1;
	  filters = 0x49494949;
	}
	break;
      case 291:
      case 50713:			/* BlackLevelRepeatDim */
	blrr = get2( fp );
	blrc = get2( fp );
	break;
      case 61450:
	blrr = blrc = 2;
      case 50714:			/* BlackLevel */
	black = getreal(fp,type);
	if (!filters || !~filters) break;
	dblack[0] = black;
	dblack[1] = (blrc == 2) ? getreal(fp,type):dblack[0];
	dblack[2] = (blrr == 2) ? getreal(fp,type):dblack[0];
	dblack[3] = (blrc == 2 && blrr == 2) ? getreal(fp,type):dblack[1];
	if (colors == 3)
	  filters |= ((filters >> 2 & 0x22222222) |
		      (filters << 2 & 0x88888888)) & filters << 1;
	FORC4 cblack[filters >> (c << 1) & 3] = dblack[c];
	black = 0;
	break;
      case 50715:			/* BlackLevelDeltaH */
      case 50716:			/* BlackLevelDeltaV */
	for (num=i=0; i < len; i++)
	  num += getreal(fp,type);
	black += num/len + 0.5;
	break;
      case 50717:			/* WhiteLevel */
	maximum = getint(fp,type);
	break;
      case 50718:			/* DefaultScale */
	pixel_aspect  = getreal(fp,type);
	pixel_aspect /= getreal(fp,type);
	break;
      case 50721:			/* ColorMatrix1 */
      case 50722:			/* ColorMatrix2 */
	FORCC for (j=0; j < 3; j++)
	  cm[c][j] = getreal(fp,type);
	use_cm = 1;
	break;
      case 50723:			/* CameraCalibration1 */
      case 50724:			/* CameraCalibration2 */
	for (i=0; i < colors; i++)
	  FORCC cc[i][c] = getreal(fp,type);
	break;
      case 50727:			/* AnalogBalance */
	FORCC ab[c] = getreal(fp,type);
	break;
      case 50728:			/* AsShotNeutral */
	FORCC asn[c] = getreal(fp,type);
	break;
      case 50729:			/* AsShotWhiteXY */
	xyz[0] = getreal(fp,type);
	xyz[1] = getreal(fp,type);
	xyz[2] = 1 - xyz[0] - xyz[1];
	FORC3 xyz[c] /= d65_white[c];
	break;
      case 50752:
	read_shorts (fp, cr2_slice, 3);
	break;
      case 50829:			/* ActiveArea */
	top_margin = getint(fp,type);
	left_margin = getint(fp,type);
	height = getint(fp,type) - top_margin;
	width = getint(fp,type) - left_margin;
	break;
      case 65026:
	if (type == 2) fp.gets(model2, 64);
    }
    fp.seek (save);
  }
  for (i=0; i < colors; i++)
    FORCC cc[i][c] *= ab[i];
  if (use_cm) {
    FORCC for (i=0; i < 3; i++)
      for (cam_xyz[c][i]=j=0; j < colors; j++)
	cam_xyz[c][i] += cc[c][j] * cm[j][i] * xyz[i];
    cam_xyz_coeff (cam_xyz);
  }
  if (asn[0]) {
    cam_mul[3] = 0;
    FORCC cam_mul[c] = 1 / asn[c];
  }
  if (!use_cm)
    FORCC pre_mul[c] /= cc[c][c];
  return 0;
}

int CLASS parse_tiff ( ASTAR_FILE_RAM &fp, int base)
{
  int doff;

  fp.seek( base );
  order = get2( fp );
  if (order != 0x4949 && order != 0x4d4d) return 0;
  get2( fp );
  while ((doff = get4( fp ))) {
    fp.seek (doff+base);
    if (parse_tiff_ifd (fp, base)) break;
  }
  return 1;
}

void CLASS apply_tiff( ASTAR_FILE_RAM &fp )
{
  int max_samp=0, raw=-1, thm=-1, i;
  struct jhead jh;

  thumb_misc = 16;
  if (thumb_offset) {
    fp.seek (thumb_offset);
    if (ljpeg_start (fp, &jh, 1)) {
      thumb_misc   = jh.bits;
      thumb_width  = jh.wide;
      thumb_height = jh.high;
    }
  }
  for (i=0; i < tiff_nifds; i++) {
    if (max_samp < tiff_ifd[i].samples)
	max_samp = tiff_ifd[i].samples;
    if (max_samp > 3) max_samp = 3;
    if ((tiff_ifd[i].comp != 6 || tiff_ifd[i].samples != 3) &&
	(tiff_ifd[i].width | tiff_ifd[i].height) < 0x10000 &&
	tiff_ifd[i].width*tiff_ifd[i].height > raw_width*raw_height) {
      raw_width     = tiff_ifd[i].width;
      raw_height    = tiff_ifd[i].height;
      tiff_bps      = tiff_ifd[i].bps;
      tiff_compress = tiff_ifd[i].comp;
      data_offset   = tiff_ifd[i].offset;
      tiff_flip     = tiff_ifd[i].flip;
      tiff_samples  = tiff_ifd[i].samples;
      raw = i;
    }
  }
  for (i=tiff_nifds; i--; )
    if (tiff_ifd[i].flip) tiff_flip = tiff_ifd[i].flip;
  if (raw >= 0 && !load_raw)
    switch (tiff_compress) {
      case 0:  case 1:
	switch (tiff_bps) {
	  case  8: assert( 0 ); break;
	  case 12: assert( 0 ); break;
	  case 14:
	  case 16: assert(0); break;
	}
	if (tiff_ifd[raw].bytes*5 == raw_width*raw_height*8) {
		assert( 0 );
	}
	break;
      case 6:  case 7:  case 99:
	load_raw = &CLASS lossless_jpeg_load_raw;		break;
      case 262:
        assert( 0 );
      case 32767:
        assert( 0 );
      case 32769:
	load_flags++;
      case 32770:
      case 32773:
        assert( 0 );
      case 34713:
        assert( 0 );
      case 65535:
        assert( 0 );
      case 65000:
	switch (tiff_ifd[raw].phint) {
		assert( 0 );
	}
      case 32867: break;
      default: is_raw = 0;
    }
  if (!dng_version)
    if ( (tiff_samples == 3 && tiff_ifd[raw].bytes &&
	  tiff_bps != 14 && tiff_bps != 2048)
      || (tiff_bps == 8 && !strstr(make,"KODAK") && !strstr(make,"Kodak") &&
	  !strstr(model2,"DEBUG RAW")))
      is_raw = 0;
  for (i=0; i < tiff_nifds; i++)
    if (i != raw && tiff_ifd[i].samples == max_samp &&
	tiff_ifd[i].width * tiff_ifd[i].height / SQR(tiff_ifd[i].bps+1) >
	      thumb_width *       thumb_height / SQR(thumb_misc+1)) {
      thumb_width  = tiff_ifd[i].width;
      thumb_height = tiff_ifd[i].height;
      thumb_offset = tiff_ifd[i].offset;
      thumb_length = tiff_ifd[i].bytes;
      thumb_misc   = tiff_ifd[i].bps;
      thm = i;
    }
  if (thm >= 0) {
    thumb_misc |= tiff_ifd[thm].samples << 5;
    switch (tiff_ifd[thm].comp) {
      case 0:
	break;
      case 1:
	break;
      case 65000:
        assert( 0 );
    }
  }
}

/*
   Identify which camera created this file, and set global variables
   accordingly.
 */
void CLASS identify( ASTAR_FILE_RAM &fp )
{
  char head[32];
  int i, c, is_canon;

  static const char *corp[] = { "Canon" };

  tiff_flip = filters = -1;	/* 0 is valid, so -1 is unknown */
  raw_height = raw_width = fuji_width = fuji_layout = cr2_slice[0] = 0;
  maximum = height = width = top_margin = left_margin = 0;
  cdesc[0] = desc[0] = artist[0] = make[0] = model[0] = model2[0] = 0;
  iso_speed = shutter = aperture = focal_len = unique_id = 0;
  tiff_nifds = 0;
  memset (tiff_ifd, 0, sizeof tiff_ifd);
  memset (gpsdata, 0, sizeof gpsdata);
  memset (cblack, 0, sizeof cblack);
  memset (white, 0, sizeof white);
  thumb_offset = thumb_length = thumb_width = thumb_height = 0;
  load_raw = 0;
  data_offset = meta_length = tiff_bps = tiff_compress = 0;
  kodak_cbpp = zero_after_ff = dng_version = load_flags = 0;
  timestamp = shot_order = tiff_samples = black = is_foveon = 0;
  mix_green = profile_length = data_error = zero_is_bad = 0;
  pixel_aspect = is_raw = raw_color = 1;
  tile_width = tile_length = INT_MAX;
  for (i=0; i < 4; i++) {
    cam_mul[i] = i == 1;
    pre_mul[i] = i < 3;
    FORC3 cmatrix[c][i] = 0;
    FORC3 rgb_cam[c][i] = c == i;
  }
  colors = 3;
  for (i=0; i < 0x4000; i++) curve[i] = i;

  order = get2( fp );
  get4( fp );
  fp.seek (0);
  fp.read(head, 32);
  fp.seek( fp.size());
  //fseek (ifp, 0, SEEK_END);
  //ftell(ifp);
  if (order == 0x4949 || order == 0x4d4d) {
    if (parse_tiff(fp, 0)) apply_tiff( fp );
  }

  for (i=0; i < sizeof corp / sizeof *corp; i++)
    if (strstr (make, corp[i]))		/* Simplify company names */
	strcpy (make, corp[i]);

  desc[511] = artist[63] = make[63] = model[63] = model2[63] = 0;

  if (!height) height = raw_height;
  if (!width)  width  = raw_width;

  if ((is_canon = !strcmp(make,"Canon")))
  {
    load_raw = &CLASS lossless_jpeg_load_raw;
  }

/* Set parameters based on camera name (for non-DNG files). */

  if (is_canon && raw_width == 5344) {
    top_margin = 51;
    left_margin = 142;
    height -= top_margin;
    width  -= left_margin;
  }

  if (!model[0])
    sprintf (model, "%dx%d", width, height);
  if (filters == UINT_MAX) filters = 0x94949494;
  if (!tiff_bps) tiff_bps = 12;
  if (!maximum) maximum = (1 << tiff_bps) - 1;
  if (!load_raw || height < 22) is_raw = 0;
  if (!cdesc[0])
    strcpy (cdesc, colors == 3 ? "RGBG":"GMCY");
  if (!raw_height) raw_height = height;
  if (!raw_width ) raw_width  = width;
  if (filters && colors == 3)
    filters |= ((filters >> 2 & 0x22222222) |
		(filters << 2 & 0x88888888)) & filters << 1;
}

void aimage_convert( AIMAGE &dst_img )
{
#ifdef GONE
	int y;

	for ( y = 0; y < iheight; ++y )
	{
		int x;
		for( x = 0; x < iwidth; ++x )
		{
			int soff = x + y * iwidth;

			unsigned int r, g, b;
			APIXEL p;

			r = image[ soff ][ 0 ];
			g = image[ soff ][ 1 ];
			b = image[ soff ][ 2 ];

			r = r * FULL_INT >> 16;
			g = g * FULL_INT >> 16;
			b = b * FULL_INT >> 16;

			p.r = r;
			p.g = g;
			p.b = b;

			dst.setpixel( x, y, p );
		}
	}
#endif
	int width, height;

	width = iwidth; 
	height = iheight; 

	assert( ( width % 2 ) == 0 );

	if ( (height % 2) == 1 )
	{
		height--;
	}

	dst_img.alloc_image( width / 2, height / 2 );

	AIMAGE::iterator dst;
	AIMAGE::iterator dst_end;

	ushort (*src_0)[4];
	ushort (*src_1)[4];

	dst = dst_img.begin();
	dst_end = dst_img.end();

	src_0 = image;
	src_1 = src_0;
	src_1 += width;
	int xp = 0;

	while ( dst != dst_end )
	{
		APIXEL p;

		p.r = (*src_0)[0];
		++src_0;
		p.g = (*src_0)[1];
		++src_0;
		++src_1;
		p.b = (*src_1)[2];
		++src_1;
		xp+=2;

		//fprintf(stderr, "%d %d %d\n", p.r, p.g, p.b );

		if ( xp == width )
		{
			src_0 += width;
			src_1 += width;
			xp = 0;
		}  
		*dst = p;
		++dst;
	}
}

void AIMAGE_DCRAW_UTILS::load(AIMAGE &dst, const char *fname )
{
  int status;
  int verbose=0;
	char cache_file[ 256 ];
	sprintf( cache_file, ".cache/%s.ppm", fname );

	FILE *f;
	if ( ( f = fopen( cache_file, "rb" )) != NULL )
	{
		fprintf(stderr,"Cache hit in raw file load\n");
		fclose( f );
		AIMAGE_PPM_UTILS::load( dst, cache_file );
		fprintf(stderr,"done load from cache\n");
		return;
	}

  ASTAR_FILE_DISK fpd( fname );
  if ( verbose ) fprintf(stderr, "Reading ram disk for raw file\n");
  ASTAR_FILE_RAM fp( fpd );
  if ( verbose ) fprintf(stderr, "Done reading ram disk for raw file\n");

#ifdef GONE
  ifname = fname;
  if (!(ifp = fopen (ifname, "rb"))) {
    perror (ifname);
    exit( 1 );
  }
#endif
  status = (identify( fp ),!is_raw);

  assert( !status );

  iheight = height;
  iwidth  = width;

  if ( verbose ) fprintf(stderr, "Alloc %d %d Meta %d DataO %d\n", iwidth , iheight, meta_length, (int) data_offset );

  image = (ushort (*)[4]) calloc (iheight*iwidth, sizeof *image);

  if (meta_length) {
  	meta_data = (char *) malloc (meta_length);
	merror (meta_data, "main()");
  }

   //fseeko (ifp, data_offset, SEEK_SET);
   fp.seek( data_offset );
  if ( verbose ) fprintf(stderr, "Loading Raw\n");
  (*load_raw)(fp);
  if ( verbose ) fprintf(stderr, "Scaling Colors\n");
  //scale_colors();
  if ( verbose ) fprintf(stderr, "Converting to Aimage\n");
  aimage_convert( dst );
  if ( verbose ) fprintf(stderr, "Done raw to aimage\n");

  free( image );

	system("mkdir .cache");
		AIMAGE_PPM_UTILS::save( dst, cache_file, true );
	
}

void AIMAGE_DCRAW_UTILS::torgb(AIMAGE &dst_img, const AIMAGE &src_img )
{
	int width, height;

	src_img.get_dim( width, height );

	assert( ( width % 2 ) == 0 );

	if ( (height % 2) == 1 )
	{
		height--;
	}

	dst_img.alloc_image( width / 2, height / 2 );

	AIMAGE::iterator dst;
	AIMAGE::iterator dst_end;
	AIMAGE::const_iterator src_0;
	AIMAGE::const_iterator src_1;

	dst = dst_img.begin();
	dst_end = dst_img.end();

	src_0 = src_img.const_begin();
	src_1 = src_0;
	src_1 += width;
	int xp = 0;

	while ( dst != dst_end )
	{
		APIXEL p;

		p.r = (*src_0).r;
		++src_0;
		p.g = (*src_0).g;
		++src_0;
		++src_1;
		p.b = (*src_1).b;
		++src_1;
		xp+=2;

		//fprintf(stderr, "%d %d %d\n", p.r, p.g, p.b );

		if ( xp == width )
		{
			src_0 += width;
			src_1 += width;
			xp = 0;
		}  
		*dst = p;
		++dst;
	}
}


